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ABSTRACT 

The  performance  of  density  functional  theories 
(DFT)  in  predicting  structural  parameters  for  six 
conventional  energetic  materials  (EM)  over  various 
degrees  of  compression  was  examined  for  a  wide  range  of 
pressures.  The  systems  studied  were  nitromethane, 
l,3,5,7-tetranitro-l,3,5,7-tetraazacyclooctane  (HMX), 
cyclotrimethylenetrinitramine  (RDX),  2,4,6,8,10,12- 
hexanitrohexaazaisowurzitane  (CL-20),  2,4,6-trinitro- 
1,3,5-benzenetriamine  (TATB),  and  pentaerythritol 
tetranitrate  (PETN).  Dependencies  of  results  on  basis  set 
size,  k-points,  choice  of  pseudopotential  and  method  were 
explored.  The  results  indicate  that  at  zero  compression, 
DFT  is  not  adequate  to  describe  crystallographic 
parameters  for  the  systems  under  study.  Flowever,  at 
compressions  consistent  with  6  GPa  or  greater,  DFT 
predictions  of  crystal  volumes  are  within  5%  of 
experiment,  and  are  insensitive  to  method,  choice  of 
pseudopotential  and  basis  set  size.  The  results  suggest 
that  the  major  source  of  error  in  DFT  calculations  applied 
to  systems  similar  to  these  are  due  to  inadequate  treatment 
of  van  der  Waals  forces,  which  are  the  dominant  forces  in 
molecular  organic  crystals  at  the  ambient  state. 

1.  INTRODUCTION 

Up  to  the  last  decade,  available  computing  processing 
power  limited  theoretical  treatment  of  condensed  phase 
systems  to  classical  molecular  simulations  methods  such 
as  molecular  dynamics  or  Monte  Carlo.  While  useful  in 
predicting  equation  of  state  (EOS)  information  or 
providing  temporal  descriptors  of  the  response  of  a 
material  to  stimuli,  the  accuracy  of  the  results  were  almost 
completely  dependent  on  the  model  used  to  describe  the 
interatomic  interactions.  Additionally,  such  simulations 
cannot  capture  quantum  effects  such  as  tunneling,  zero 
point  energy,  and  electronic  excitation.  On  the  other 
hand,  quantum  mechanical  methods  do  not  suffer  from 
these  limitations,  yet  have  been  precluded  from  wide¬ 
spread  use,  particularly  in  condensed  phase  systems,  due 
to  their  often  inordinate  computational  requirements.  In 
fact,  even  with  the  processing  power  available  today, 
quantum  mechanical  treatments  of  condensed  phase 
systems  are  largely  restricted  to  the  quantum  mechanical 
method  known  as  Density  Functional  Theory  (DFT). 
DFT  is  routinely  used  to  characterize  large  polyatomic 
molecules  and  certain  condensed  phase,  materials  (e.g. 


metals,  covalently  bound  systems).  For  those  cases,  DFT 
predictions  demonstrate  a  high  degree  of  accuracy. 
Unfortunately,  DFT  has  also  demonstrated  some 
spectacular  failures  in  predictive  capabilities  due  to 
approximations  in  the  implementation  of  the  Density 
Functional  Theory  [Byrd  et  al.,  2004].  Thus,  evaluation 
of  the  suitability  of  this  methodology  when  applied  to 
material  types  for  which  extensive  assessment  has  not 
been  performed  is  required.  Conventional  CHNO 
energetic  materials,  which  are  organic  molecular  crystals, 
fall  within  this  category  of  materials  for  which  assessment 
of  DFT  applications  is  necessary.  We  present  such  an 
assessment  through  a  review  of  our  work  which  explores 
the  suitability  of  various  DFT  methods  applied  to  CHNO 
molecular  crystals  and  their  performance  at  different 
degrees  of  compression.  [Byrd  and  Rice,  2007] 

2.  COMPUTATIONAL  DETAILS 

This  work  explores  basis  set,  pseudopotential, 
method  and  k-point  dependencies  on  predictions  of 
crystal  and  molecular  structural  parameters  for  six  well- 
studied  conventional  CHNO  explosives  under  different 
degrees  of  compression.  The  systems  are  RDX,  HMX, 
CL-20,  nitromethane,  PETN,  and  TATB.  The  GGA  DFT 
Perdew-Burke-Ernzerhof  (PBE)  [Perdew  et  al.,  1996] 
functional  was  also  applied  to  nitromethane  and  RDX. 
All  were  subjected  to  the  generalized  gradient 
approximation  density  functional  theory  (GGA  DFT) 
Perdew- Wang  91  (PW91)  [Perdew,  1991]  functional  as 
implemented  in  the  Vienna  Ab-Initio  Package  (VASP) 
[Kresse  and  Furthermuller,  2003].  The  Vanderbilt 
ultrasoft  pseudopotentials  (USP)  [Vanderbilt,  1990]  and 
Monkhorst-Pack  k-point  generation  method  were  used  for 
all  calculations.  Plane  wave  kinetic  energy  cut-offs  (Ecut) 
ranged  from  280  eV  to  800  eV.  Initial  crystal  structures 
used  in  the  geometry  optimizations  correspond  to 
experimental  structures.  In  the  geometry  optimizations 
space  group  symmetry  was  not  imposed  and  all  ionic  and 
cell  degrees  of  freedom  were  allowed  to  vary.  The 
electronic  energies  were  converged  to  2.0xl0'6  eV,  while 
the  structures  were  converged  once  the  difference  in  free 
energy  between  gradient  steps  was  less  than  2.0xl0"5  eV. 
Smearing  width  was  set  to  0.0001  in  the  Methfessel- 
Paxton  [Methfessel  and  Paxton,  1989]  method  minimized 
electron  smearing,  since  molecular  crystals  are  insulators. 
Also,  the  RMM-DIIS  method  was  employed  to  accelerate 
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convergence  except  for  those  cases  in  which  conjugate 
gradients  were  required  to  achieve  convergence. 

At  zero  pressure,  basis  set  and  k-point  dependencies 
were  explored  for  only  four  of  the  six  systems 
(nitromethane,  HMX,  RDX  and  CL-20).  At  higher 
degrees  of  compression,  planewave  kinetic  energy  cut¬ 
offs  (Ecut)  were  restricted  to  either  396  eV  and/or  545  eV. 
Also  at  the  higher  degrees  of  compression,  the  Local 
Density  Approximation  (LDA)  [Ceperley  and  Alder, 
1980]  functional  was  applied  to  TATB,  PETN,  and  RDX. 
Because  of  concern  over  the  adequacy  of  ultrasoft 
pseudopotentials  for  crystals  under  pressure,  we  used 
projected  augmented  wavefunction  (PAW)  [Blochl,  1994; 
Kresse  and  Joubert,  1999]  pseudopotentials  to  check  for 
possible  pseudopotential  induced  errors  for  the  TATB  and 
PETN  crystals.  To  further  eliminate  pseudopotentials 
from  consideration  as  the  root  of  any  error,  norm- 
conserving  Troullier-Martins  pseudopotentials  and  the 
PBE  functional  as  implemented  in  the  PWSCF  program 
package  [Baroni  et  al.,  2008]  were  used  to  predict  crystal 
structures  of  compressed  PETN.  For  these  calculations, 
optimized  geometries  were  obtained  using  a  damped 
molecular  dynamics  calculations,  with  convergence 
criteria  for  electron  energies  set  at  1.0x1  O'9  Rydberg  and 
total  energy  convergence  set  at  1.0x1  O'6  Rydberg.  The 
final  converged  geometries  were  obtained  through  a 
damped  molecular  dynamics  calculation  as  coded  in 
PWSCF.  k-point  meshes  employed  for  non-zero 
pressures  were:  2x1x2  for  HMX,  lxlxl  for  RDX  and  CL- 
20,  2x2x2  for  TATB,  and  2x2x2  for  PETN.  The 
Troullier-Martins  PBE  PWSCF  calculations  used  only  the 
gamma  point.  Initial  structures  used  in  the  geometry 
optimizations  of  all  crystals  under  all  degrees  of 
compression  correspond  to  the  ambient  pressure 
experimental  crystal  structure. 


3.  RESULTS  AND  DISCUSSION 
3.1  Predictions  for  crystals  with  no  compression 

Starting  with  investigating  the  zero  compression 
crystals,  for  each  of  the  four  crystals  for  which  basis  set 
and  k-point  dependencies  were  explored  (nitromethane 
[Trevino  et  al.,  1980],  FIMX  [Choi  and  Boutin,  1970], 
RDX  [Choi  and  Prince,  1972],  CL-20  [Nielsen  et  al, 
1998]),  the  grid  that  had  the  minimum  number  of  k-points 
that  gave  a  converged  result  was  identified.  Once  the  grid 
was  identified,  crystallographic  and  molecular  structural 
parameters  were  calculated  for  kinetic  energy  cutoffs 
ranging  between  280  and  800  eV.  The  resulting  grids  for 
nitromethane  and  HMX  are  (2x2x2)  and  (2x1x2), 
respectively.  For  RDX  and  CL-20,  we  determined  the 
(lxlxl)  was  sufficient  to  yield  converged  results. 


Predicted  molecular  structural  parameters  (bond 
lengths  and  bond  angles)  were  converged  between  495 
and  545  eV,  and  once  converged,  are  in  excellent 
agreement  with  measured  values.  Mean  errors  in  bond 
lengths  for  all  four  systems  are  no  more  than  0.017  A, 
with  the  mean  error  in  bond  lengths  for  CL-20  an  order  of 
magnitude  smaller.  The  maximum  errors  using  a 
sufficient  basis  set  for  bond  distances  is  less  than  0.04  A 
for  nitromethane,  RDX  and  HMX,  and  0.02  A  for  CL-20. 
Similar  to  predictions  for  isolated  molecules,  the  DFT 
predictions  of  molecular  structural  parameters  are  in  very 
good  agreement  with  experimental  values.  However,  for 
the  crystallographic  structural  parameters,  we  observe 
poor  convergence  to  erroneous  values.  This  is  illustrated 
in  Figure  1,  which  shows  the  percent  error  in  lattice 
parameters  as  a  function  of  kinetic  energy  cutoff  for 
nitromethane,  HMX,  RDX  and  CL-20. 


Kinetic  Energy  Cutoff  (eV) 

Figure  1.  Variation  of  lattice  dimensions  with  pressure  for  (a) 
nitromethane,  (b)  HMX,  (c)  RDX,  and  (d)  CL-20.  The  a,  b,  C- 
cell  constants  are  represented  by  circles,  open  circles  and 
triangles,  respectively. 

The  range  of  kinetic  energy  cutoffs  for  three  of  the 
systems  is  280  to  700  eV;  an  additional  point  was 
calculated  for  RDX  at  800  eV  since  it  was  not  clear  that 
the  results  had  converged  at  700  eV.  It  is  interesting  to 
note  that  at  around  400  eV,  near  the  default  cutoff  value 
in  VASP  for  an  oxygen  containing  compound  (395  eV), 
the  percent  error  in  lattice  parameters  for  all  systems  is 
approximately  zero.  However,  the  results  do  not  begin  to 
converge  until  545  eV,  but  do  so  to  values  between  5  to 


2 


10%  too  large.  Although  not  shown  in  these  figures,  we 
observe  similar  behavior  for  TATB  and  PETN.  We 
should  point  out  that  space  group  symmetry  is  conserved 
in  the  geometry  optimization,  and  for  the  two  non- 
orthorhombic  cells  (CL-20  and  HMX),  the  deviation  of 
the  non-90°  angle  is  within  1.5°  of  experiment. 

These  errors  in  the  lattice  vectors  could  be  inherent  to 
the  PW91  functional.  To  explore  this  possibility,  we 
calculated  crystallographic  parameters  using  the  PBE 
functional  as  implemented  in  the  VASP  suite  using  Ecut  of 
495,  545,  and  645  eV  for  nitromethane  and  Ecut  of  430 
and  495  eV  for  RDX.  Figure  2  illustrates  the  errors  in  the 
a,  b,  and  c  lattice  vectors  for  both  nitromethane  and  RDX 
with  this  method;  results  using  the  PW91  functional  are 
provided  for  comparison.  As  evident  in  the  figure, 
although  the  results  using  the  PBE  functional  are  smaller 
than  those  generated  using  the  PW91  functional,  the 
improvement  is  not  substantial  enough  to  conclude  that 
the  major  source  of  error  is  in  the  choice  of  the  PW91 
functional.  More  likely,  the  major  source  of  error  is  due 
to  the  inherent  deficiencies  in  DFT  to  appropriately  model 
van  der  Waals  forces,  which  are  the  main  binding  forces 
in  molecular  organic  crystals  such  as  these. 


Cell  Edge  Length 


Cell  Edge  Length 


Figure  2.  Comparison  of  percent  errors  versus  experiment  for 
cell  lengths  between  PW91  and  PBE  for  (a)  nitromethane  and 
(b)  RDX 

3.2  Predictions  at  varying  degrees  of  compression 

Although  DFT  is  clearly  unable  to  adequately  predict 
crystal  structures  of  energetic  molecular  crystals  at  zero 


compression,  there  are  studies  indicating  weakly-bound 
molecular  systems  under  sufficient  compression  can  be 
reasonably  described  by  DFT  [Miao  et  al.,  2003].  This, 
as  well  as  our  interests  in  theoretical  characterization  of 
energetic  materials  under  extreme  conditions  (including 
high  compression),  we  explored  the  ability  of 
conventional  DFT  to  predict  structures  of  energetic 
materials  under  high  compression.  Also,  we  were 
interested  in  determining  the  degree  of  compression  for 
which  DFT  could  suitably  treat  the  EM;  i.e.  those 
pressure  regimes  for  which  dispersion  is  not  the  dominant 
component  of  the  interactions.  Therefore,  we  calculated 
crystallographic  and  intramolecular  parameters  for  five 
energetic  molecular  crystals  under  varying  degrees  of 
compression.  In  this  paper,  we  report  results  only  for  the 
systems  for  which  there  are  extensive  experimental 
information  over  a  large  pressure  range  (TATB  [Olinger 
and  Cady,  1974],  HMX  [Olinger  et  al.,  1978],  and  PETN 
[Olinger  et  al.,  1975]).  We  have  also  calculated  structural 
parameters  of  compressed  RDX  [Olinger  et  al.,  1978]  and 
CL-20  [Sorescu  et  al.,  1998];  however,  both  systems 
undergo  polymorphic  phase  changes  at  relatively  low 
pressures,  and  therefore,  there  are  limited  experimental 
data  available  for  comparison. 

As  in  the  study  described  in  Section  3.1,  we  explored 
dependencies  of  the  results  on  basis  set  size,  choice  of 
pseudopotential,  and  specific  implementations  of  DFT. 
With  regard  to  the  latter  point,  we  performed  calculations 
using  the  545  eV  planewave  basis  set  and  Vanderbilt 
ultrasoft  pseudopotentials  (USP)  using  LDA  and  GGA 
(PW91  and  PBE)  methods.  LDA  volumes  for  the  TATB 
and  PETN  crystals  at  zero  pressures  are  smaller  than 
experiment  by  11%  and  15%,  respectively.  As  pressures 
increase,  LDA  errors  for  TATB  and  PETN  decrease  to 
4%  too  small  at  7  GPa  for  TATB  and  7%  too  small  at  10 
GPA  for  PETN.  Unlike  the  LDA  results,  the  GGA 
methods  predict  lattice  parameters  that  are  substantially 
larger  than  the  experimental  values  at  zero  pressure,  with 
deviations  of  the  predictions  from  experiment  decreasing 
with  increasing  pressure.  This  behavior  is  apparent  in 
Fig.  3,  which  shows  the  behavior  of  the  volume  and 
lattice  parameters  of  PETN  with  increase  in  pressure  for 
the  various  methods.  For  all  calculations  at  all  pressures, 
the  space  group  symmetry  was  maintained  in  the 
completely  unconstrained  geometry  optimizations. 

Extensive  structural  analyses  of  the  contents  of  the 
unit  cells  for  all  crystals  subjected  to  compression  were  to 
determine  any  pressure-dependent  rotational  and 
translational  disorder  of  the  molecules  within  the  unit  cell. 
Because  detailed  experimental  data  about  atomic 
positions  in  the  crystals  are  lacking  for  pressures  above  1 
atm.,  we  are  limited  to  structural  comparisons  with 
experimental  crystal  structures  at  room  conditions.  Unit 
cell  degrees  of  freedom  that  will  be  compared  are  the  six 
cell  constants  (a,  b,  c,  a,  P,  y)  and  the  orientational  and 


3 


translational  degrees  of  freedom  of  each  of  the  molecules 
in  the  unit  cells.  The  molecular  translational  degrees  of 
freedom  will  be  given  in  fractional  coordinates  (sx,  sy,  sz) 
of  the  molecular  mass  centers.  Euler  angles  for  rigid- 
body  rotation  about  inertial  axes  (9, 4>,  \|/)  for  all 
molecules  in  the  unit  cell  are  used  to  describe  the 
molecular  orientational  degrees  of  freedom.  Symmetry 
constraints  were  not  imposed  in  the  geometry 
optimizations,  therefore,  the  molecules  within  the  unit  cell 
are  not  required  to  be  symmetry  equivalents.  Thus, 
molecular  translational  and  orientational  parameters  for 
each  molecule  within  the  unit  cell  were  examined. 


Pressure  (GPa) 


Pressure  (GPa) 


Figure  3.  Variation  of  volume  and  lattice  dimensions  with 
pressure  for  PETN. 

For  RDX,  fractional  coordinates  of  the  molecular 
mass  centers  differ  slightly  from  those  in  the  ambient- 
state  crystal,  with  differences  in  the  thousandths. 
Orientation  of  the  molecules  about  any  one  of  the  Euler 
angles  is  no  greater  than  3°,  resulting  in  a  root-mean 
square  (rms)  deviation  of  the  predicted  atomic 
displacements  from  the  experimental  positions  of  0.064  A 
and  a  maximum  atomic  displacement  is  0.229  A. 
Molecular  structural  parameters  for  e-CL-20  generated  at 
2.5  GPa  are  in  similar  close  agreement  with  the  ambient 
state  crystal. 

The  differences  in  fractional  coordinates  of  the 
molecular  mass  centers  for  (3-FIMX  and  TATB  are  similar 
to  those  of  RDX  up  to  4  GPa,  but  deviations  of  the 
fractional  coordinates  become  as  large  as  0.011  at  higher 
degrees  of  compression.  Deviations  in  the  Euler  angles 
are  no  greater  than  2°  over  the  pressure  ranges.  For 
TATB,  the  rms  deviations  in  atomic  displacements  are  all 
~  0.1  A,  and  are  attributed  to  the  difference  in  the 
orientation  of  the  amine  groups  relative  to  the  ring.  In 
experiment,  the  amine  hydrogens  are  out  of  the  plane  of 
the  molecule,  whereas  in  the  theoretical  calculations,  all 
atoms  lie  in  approximately  within  a  plane. 


Four  of  the  five  systems  studied  are  relatively  rigid 
molecules,  three  of  which  are  cyclic  nitramines  and  one  a 
nitroaromatic.  The  remaining  molecule,  PETN,  however, 
is  an  acyclic,  branched-chain  molecule  that  would  more 
readily  deform  upon  compression.  This  expectation  is 
supported  by  results  from  diamond  anvil  cell  (DAC) 
experiments  indicating  a  conformation  change  to  a  lower- 
symmetry  molecular  structure  after  subjecting  the  system 
to  pressures  greater  than  5  GPa  [Grudzkov  et  ah,  2004]. 
A  second  group  performed  a  series  of  DAC  experiments 
in  which  changes  in  PETN  with  high  pressure  were 
monitored  [Lipinska-Kalita  et  ah,  2005].  These 
experiments  were  performed  without  and  with  a  pressure- 
transmitting  medium  (nitrogen),  with  the  latter  designed 
to  provide  quasi-hydrostatic  conditions.  A  structural 
phase  transition  was  apparent  at  8  GPa  in  the  experiments 
in  which  the  pressure-transmitting  medium  was  not 
present,  but  results  in  which  nitrogen  was  introduced  into 
the  cell  indicated  that  any  structural  rearrangements  due 
to  compression  under  these  experimental  conditions  are 
fully  reversible.  These  authors  conclude  that  the 
discrepancies  between  the  two  sets  of  measurements  are 
due  to  differences  in  the  experimental  condition 
specifically  that  more  severe  non-hydrostatic  conditions 
are  introduced  in  the  absence  of  the  pressure-transmitting 
medium.  Our  calculations  represent  pure  hydrostatic 
compression;  differences  in  fractional  coordinates  are  less 
than  0.03  over  the  entire  range  of  pressures  studied,  with 
differences  in  the  Euler  angles  are  no  greater  than  3.5°. 
At  the  largest  pressure  (10.45  GPa),  the  rms  deviation  of 
the  atomic  displacement  is  0.117  A  and  the  maximum 
deviation  is  0.170A.  A  graphical  depiction  of  the 
predicted  unit  cell  at  10.45  GPa  superimposed  on  the 
experimental  cell  at  ambient  conditions  is  given  in  Fig.  4, 
demonstrating  the  small  deformation  of  the  molecular 
structure  at  high  degrees  of  compression.  The  results 
presented  here  are  supportive  of  the  suggestion  that  under 
quasi-hydrostatic  compression,  no  irreversible  structural 
changes  occur.  Flowever,  additional  experimentation  or 
theoretical  investigation  of  non-hydrostatic  compression 
needs  to  be  performed  to  provide  a  definitive  conclusion. 


Figure  4.  Pictorial  view  of  the  10.45  GPa  PW91  (545 
eV)  unit  cell  structure  of  PETN  superimposed  onto  the 
ambient  pressure  experimental  unit  cell. 
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Basis  set  size  dependencies  for  compressed  crystals 
were  also  investigated  using  the  PW91  functional  and  396 
eV  kinetic  energy  cutoff.  For  pressures  greater  than  6 
GPa,  the  predictions  using  the  396  eV  and  545  eV  basis 
sets  converge  to  the  same  values,  and  are  all  slightly 
larger  than  experiment.  Flowever,  the  error  is 
substantially  less  than  that  at  545  eV  and  zero  pressure.  A 
larger  kinetic  energy  cutoff  (645  eV)  was  utilized  in 
calculations  for  two  pressures  for  each  of  TATB  and 
PETN  to  ensure  that  the  545  eV  kinetic  energy  cutoff 
produced  converged  results;  the  difference  between  the 
predicted  volumes  at  the  various  pressures  were  less  than 
2  A3  for  all  cases. 

The  results  also  appeared  to  be  insensitive  to  choice 
of  functional.  For  pressures  greater  than  2  GPa,  PBE 
using  PAW  pseudopotentials  and  PW91  using  USP 
converge  to  the  same  result,  again  illustrated  for  PETN  in 
Figure  3.  Finally,  to  ensure  that  the  softer  PAW  and 
ultrasoft  pseudopotentials  were  not  contributing  to  error  at 
compression,  we  employed  an  extremely  hard  Troullier- 
Martins  (TM)  pseudopotential  with  a  series  of  planewave 
basis  sets  for  PETN.  The  results  indicate  that  the 
converged  TM  results  (using  extremely  large  basis  sets) 
are  comparable  to  those  generated  using  the  softer 
pseudopotentials.  As  an  added  bonus,  these  calculations, 
which  were  performed  using  the  software  suite  PWSCF 
[Baroni  et  al,  2008],  have  allowed  us  to  eliminate  as 
possible  sources  of  error  the  specific  implementations  of 
quantum  mechanical  methods  in  different  software 
packages. 


4.  FUTURE  DIRECTIONS 

We  are  in  the  process  of  evaluating  emerging 
methods  designed  to  remedy  the  deficiency  of 
conventional  DFT,  i.e.  its  improper  treatment  of 
dispersion.  To  date,  we  have  explored  a  method 
developed  by  Langreth  (vdW-DFT)  [Rydberg  et  al.,  2003; 
Dion  et  al.,  2004]  to  correct  for  the  inadequate  treatment 
of  van  der  Waals  in  DFT  [Hooper  et  al,  2008];  however, 
while  yielding  improved  results  over  conventional  DFT, 
the  approach  was  computational  unwieldy.  We  also 
exploring  other  alternatives,  including  new  functionals 
developed  by  University  of  Minnesota  (the  M0X  series) 
[Zhao  et  al.,  2005;  Zhao  et  al.,  2006;  Zhao  and  Truhlar, 
2008a,  2008b]  and  Dispersion  Corrected  Atom  Centered 
Pseudopotentials  (DCACPs)  [von  Lilienfeld  et  al.,  2004; 
von  Lilienfeld  et  al.,  2005;  Tkatchenko  and  von 
Lilienfeld,  2006;  von  Lilienfeld  and  Andrienko,  2006]. 

Initially,  the  M0X  functions  were  not  employed  in 
any  periodic  DFT  codes,  but  only  in  molecular  codes; 
thus,  we  assessed  the  functionals  for  various  dimers  while 
concurrently  augmenting  CP2K  (see  discussion  on  CP2K 
above)  to  accommodate  the  M0X  family  of  functionals. 


While  efforts  are  ongoing  to  implement  M0X  functionals 
into  CP2K,  initial  results  for  this  family  of  functionals 
[Hooper,  2008]  suggest  these  functions  are  promising. 
Figure  5  shows  the  binding  energy  of  an  RDX  dimer 
using  the  M05  functional,  the  Langreth  method  (vdW- 
DFT),  conventional  PBE  DFT  and  the  highly  accurate 
(and  extremely  computationally  expensive)  SAPT 
method.  Not  surprisingly,  the  standard  PBE  DFT 
performs  poorly  when  compared  to  SAPT,  while  the  other 
DFT  methods  are  more  accurate.  While  the  vdW-DFT  is 
in  closer  agreement  with  the  SAPT  value,  one  must  take 
into  account  that  this  method  requires  30-50%  more 
calculation  time  per  point  than  the  M05  method.  As  the 
more  advanced  M06  models  are  coded,  we  expect  to  see 
an  improvement  over  the  results  predicted  with  M05. 


Radius  (angstrom) 

Figure  5.  Binding  Energy  as  a  function  of  RDX  dimer 
separation  [Hooper,  2008]. 

We  are  in  the  process  of  investigating  the  accuracy  of 
DFT  calculations  of  molecular  crystals  using  DCACPs,  a 
capability  developed  to  include  van  der  Waals 
interactions  in  a  self-consistent  manner  within  the 
framework  of  DFT.  This  capability  was  present  in  CP2K, 
and  we  have  begun  an  exploration  into  its  performance 
when  applied  to  HMX,  RDX,  TATB  and  PETN  at  zero 
compression  (where  conventional  DFT  fails)  as  well  as  at 
the  higher-pressure  regimes  where  conventional  DFT 
better  performs.  The  results  will  be  compared  both  with 
experimental  information  and  with  calculations  using 
conventional  DFT. 

Calculations  to  date  have  only  explored  the 
dependence  of  the  results  on  basis  set  size.  All 
calculations  were  performed  using  TZVP  basis  sets  and 
the  PBE  functional  using  DCACPs.  We  intend  to  perform 
a  series  of  calculations  for  those  systems  for  which 
conventional  DFT  fails,  i.e.  HMX,  PETN,  TATB  and 
RDX.  Preliminary  results  for  PETN  are  shown  in  Fig.  6. 
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Additionally  calculations  for  the  same  system  using 
conventional  DFT  and  the  PBE  functional  are  given  for 
comparison.  it  is  clear  that  for  PETN,  the  PBE 
functional  using  DCACPs  performs  quite  well  in 
predicting  the  low  pressure  crystal  structure. 


Figure  6.  Unit  Cell  volumes  versus  plane  wave  cutoff  for 
PETN.  Black  circles  and  lines  denote  conventional  DFT 
results;  red  circles  and  lines  denote  DCACP  DFT  results 


CONCLUSIONS 

We  have  presented  a  review  of  our  efforts  to  evaluate 
the  performance  of  various  density  functional  theories 
when  applied  to  a  few  common  energetic  materials  both 
at  zero  pressure  and  under  compression.  The  results  show 
that  DFT  is  not  suitable  for  accurate  prediction  of 
crystallographic  parameters  of  these  systems  at  low 
degrees  of  compression;  percent  differences  from 
experimental  volumes  range  from  10  to  16%  for  the 
systems  studied.  However,  molecular  structural 
parameters  at  zero  pressure  are  in  good  agreement  with 
experimental  values.  At  higher  pressures,  PW91  and  PBE 
functionals  converged  to  the  same  results  regardless  of 
basis  set  or  pseudopotential  and  over-predict  experimental 
volumes  by  approximately  5%. 

These  results  support  the  hypothesis  that  the  majority 
of  the  error  in  these  calculation  at  low  compression  due  to 
inadequate  treatment  of  the  weak  binding  forces  in 
molecular  crystal  (i.e.  van  der  Waals  forces).  However, 
the  results  clearly  indicate  that  for  compression  regimes  in 
which  other  intermolecular  forces  dominate,  DFT  is  a 
suitable  method  to  predict  crystallographic  parameters. 
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